
<!DOCTYPE html>

<html lang="en">
  <head>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0" />
    <title>tigramite.toymodels.structural_causal_processes &#8212; Tigramite 5.2 documentation</title>
    <link rel="stylesheet" type="text/css" href="../../../_static/pygments.css" />
    <link rel="stylesheet" type="text/css" href="../../../_static/alabaster.css" />
    <script data-url_root="../../../" id="documentation_options" src="../../../_static/documentation_options.js"></script>
    <script src="../../../_static/jquery.js"></script>
    <script src="../../../_static/underscore.js"></script>
    <script src="../../../_static/_sphinx_javascript_frameworks_compat.js"></script>
    <script src="../../../_static/doctools.js"></script>
    <link rel="index" title="Index" href="../../../genindex.html" />
    <link rel="search" title="Search" href="../../../search.html" />
   
  <link rel="stylesheet" href="../../../_static/custom.css" type="text/css" />
  
  
  <meta name="viewport" content="width=device-width, initial-scale=0.9, maximum-scale=0.9" />

  </head><body>
  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          

          <div class="body" role="main">
            
  <h1>Source code for tigramite.toymodels.structural_causal_processes</h1><div class="highlight"><pre>
<span></span><span class="sd">&quot;&quot;&quot;Tigramite toymodels.&quot;&quot;&quot;</span>

<span class="c1"># Author: Jakob Runge &lt;jakob@jakob-runge.com&gt;</span>
<span class="c1">#</span>
<span class="c1"># License: GNU General Public License v3.0</span>
<span class="kn">from</span> <span class="nn">__future__</span> <span class="kn">import</span> <span class="n">print_function</span>
<span class="kn">from</span> <span class="nn">collections</span> <span class="kn">import</span> <span class="n">defaultdict</span><span class="p">,</span> <span class="n">OrderedDict</span>
<span class="kn">import</span> <span class="nn">sys</span>
<span class="kn">import</span> <span class="nn">warnings</span>
<span class="kn">import</span> <span class="nn">copy</span>
<span class="kn">import</span> <span class="nn">math</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">scipy.sparse</span>
<span class="kn">import</span> <span class="nn">scipy.sparse.linalg</span>
<span class="kn">from</span> <span class="nn">numba</span> <span class="kn">import</span> <span class="n">jit</span>
<span class="kn">import</span> <span class="nn">itertools</span>

<span class="k">def</span> <span class="nf">_generate_noise</span><span class="p">(</span><span class="n">covar_matrix</span><span class="p">,</span> <span class="n">time</span><span class="o">=</span><span class="mi">1000</span><span class="p">,</span> <span class="n">use_inverse</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Generate a multivariate normal distribution using correlated innovations.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    covar_matrix : array</span>
<span class="sd">        Covariance matrix of the random variables</span>
<span class="sd">    time : int</span>
<span class="sd">        Sample size</span>
<span class="sd">    use_inverse : bool, optional</span>
<span class="sd">        Negate the off-diagonal elements and invert the covariance matrix</span>
<span class="sd">        before use</span>

<span class="sd">    return_eigenvectors</span>
<span class="sd">    -------</span>
<span class="sd">    noise : array</span>
<span class="sd">        Random noise generated according to covar_matrix</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Pull out the number of nodes from the shape of the covar_matrix</span>
    <span class="n">n_nodes</span> <span class="o">=</span> <span class="n">covar_matrix</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="c1"># Make a deep copy for use in the inverse case</span>
    <span class="n">this_covar</span> <span class="o">=</span> <span class="n">covar_matrix</span>
    <span class="c1"># Take the negative inverse if needed</span>
    <span class="k">if</span> <span class="n">use_inverse</span><span class="p">:</span>
        <span class="n">this_covar</span> <span class="o">=</span> <span class="n">copy</span><span class="o">.</span><span class="n">deepcopy</span><span class="p">(</span><span class="n">covar_matrix</span><span class="p">)</span>
        <span class="n">this_covar</span> <span class="o">*=</span> <span class="o">-</span><span class="mi">1</span>
        <span class="n">this_covar</span><span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">diag_indices_from</span><span class="p">(</span><span class="n">this_covar</span><span class="p">)]</span> <span class="o">*=</span> <span class="o">-</span><span class="mi">1</span>
        <span class="n">this_covar</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">this_covar</span><span class="p">)</span>
    <span class="c1"># Return the noise distribution</span>
    <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">multivariate_normal</span><span class="p">(</span><span class="n">mean</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n_nodes</span><span class="p">),</span>
                                            <span class="n">cov</span><span class="o">=</span><span class="n">this_covar</span><span class="p">,</span>
                                            <span class="n">size</span><span class="o">=</span><span class="n">time</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">_check_stability</span><span class="p">(</span><span class="n">graph</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Raises an AssertionError if the input graph corresponds to a non-stationary</span>
<span class="sd">    process.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    graph : array</span>
<span class="sd">        Lagged connectivity matrices. Shape is (n_nodes, n_nodes, max_delay+1)</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the shape from the input graph</span>
    <span class="n">n_nodes</span><span class="p">,</span> <span class="n">_</span><span class="p">,</span> <span class="n">period</span> <span class="o">=</span> <span class="n">graph</span><span class="o">.</span><span class="n">shape</span>
    <span class="c1"># Set the top section as the horizontally stacked matrix of</span>
    <span class="c1"># shape (n_nodes, n_nodes * period)</span>
    <span class="n">stability_matrix</span> <span class="o">=</span> \
        <span class="n">scipy</span><span class="o">.</span><span class="n">sparse</span><span class="o">.</span><span class="n">hstack</span><span class="p">([</span><span class="n">scipy</span><span class="o">.</span><span class="n">sparse</span><span class="o">.</span><span class="n">lil_matrix</span><span class="p">(</span><span class="n">graph</span><span class="p">[:,</span> <span class="p">:,</span> <span class="n">t_slice</span><span class="p">])</span>
                             <span class="k">for</span> <span class="n">t_slice</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">period</span><span class="p">)])</span>
    <span class="c1"># Extend an identity matrix of shape</span>
    <span class="c1"># (n_nodes * (period - 1), n_nodes * (period - 1)) to shape</span>
    <span class="c1"># (n_nodes * (period - 1), n_nodes * period) and stack the top section on</span>
    <span class="c1"># top to make the stability matrix of shape</span>
    <span class="c1"># (n_nodes * period, n_nodes * period)</span>
    <span class="n">stability_matrix</span> <span class="o">=</span> \
        <span class="n">scipy</span><span class="o">.</span><span class="n">sparse</span><span class="o">.</span><span class="n">vstack</span><span class="p">([</span><span class="n">stability_matrix</span><span class="p">,</span>
                             <span class="n">scipy</span><span class="o">.</span><span class="n">sparse</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="n">n_nodes</span> <span class="o">*</span> <span class="p">(</span><span class="n">period</span> <span class="o">-</span> <span class="mi">1</span><span class="p">),</span>
                                              <span class="n">n_nodes</span> <span class="o">*</span> <span class="n">period</span><span class="p">)])</span>
    <span class="c1"># Check the number of dimensions to see if we can afford to use a dense</span>
    <span class="c1"># matrix</span>
    <span class="n">n_eigs</span> <span class="o">=</span> <span class="n">stability_matrix</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="k">if</span> <span class="n">n_eigs</span> <span class="o">&lt;=</span> <span class="mi">25</span><span class="p">:</span>
        <span class="c1"># If it is relatively low in dimensionality, use a dense array</span>
        <span class="n">stability_matrix</span> <span class="o">=</span> <span class="n">stability_matrix</span><span class="o">.</span><span class="n">todense</span><span class="p">()</span>
        <span class="n">eigen_values</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">scipy</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">eig</span><span class="p">(</span><span class="n">stability_matrix</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># If it is a large dimensionality, convert to a compressed row sorted</span>
        <span class="c1"># matrix, as it may be easier for the linear algebra package</span>
        <span class="n">stability_matrix</span> <span class="o">=</span> <span class="n">stability_matrix</span><span class="o">.</span><span class="n">tocsr</span><span class="p">()</span>
        <span class="c1"># Get the eigen values of the stability matrix</span>
        <span class="n">eigen_values</span> <span class="o">=</span> <span class="n">scipy</span><span class="o">.</span><span class="n">sparse</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">eigs</span><span class="p">(</span><span class="n">stability_matrix</span><span class="p">,</span>
                                                <span class="n">k</span><span class="o">=</span><span class="p">(</span><span class="n">n_eigs</span> <span class="o">-</span> <span class="mi">2</span><span class="p">),</span>
                                                <span class="n">return_eigenvectors</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
    <span class="c1"># Ensure they all have less than one magnitude</span>
    <span class="k">assert</span> <span class="n">np</span><span class="o">.</span><span class="n">all</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">eigen_values</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1.</span><span class="p">),</span> \
        <span class="s2">&quot;Values given by time lagged connectivity matrix corresponds to a &quot;</span><span class="o">+</span>\
        <span class="s2">&quot; non-stationary process!&quot;</span>

<span class="k">def</span> <span class="nf">_check_initial_values</span><span class="p">(</span><span class="n">initial_values</span><span class="p">,</span> <span class="n">shape</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Raises a AssertionError if the input initial values:</span>
<span class="sd">        * Are not a numpy array OR</span>
<span class="sd">        * Do not have the shape (n_nodes, max_delay+1)</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    graph : array</span>
<span class="sd">        Lagged connectivity matrices. Shape is (n_nodes, n_nodes, max_delay+1)</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Ensure it is a numpy array</span>
    <span class="k">assert</span> <span class="nb">isinstance</span><span class="p">(</span><span class="n">initial_values</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">ndarray</span><span class="p">),</span>\
        <span class="s2">&quot;User must provide initial_values as a numpy.ndarray&quot;</span>
    <span class="c1"># Check the shape is correct</span>
    <span class="k">assert</span> <span class="n">initial_values</span><span class="o">.</span><span class="n">shape</span> <span class="o">==</span> <span class="n">shape</span><span class="p">,</span>\
        <span class="s2">&quot;Initial values must be of shape (n_nodes, max_delay+1)&quot;</span><span class="o">+</span>\
        <span class="s2">&quot;</span><span class="se">\n</span><span class="s2"> current shape : &quot;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">initial_values</span><span class="o">.</span><span class="n">shape</span><span class="p">)</span><span class="o">+</span>\
        <span class="s2">&quot;</span><span class="se">\n</span><span class="s2"> desired shape : &quot;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">shape</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">_var_network</span><span class="p">(</span><span class="n">graph</span><span class="p">,</span>
                 <span class="n">add_noise</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span>
                 <span class="n">inno_cov</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span>
                 <span class="n">invert_inno</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span>
                 <span class="n">T</span><span class="o">=</span><span class="mi">100</span><span class="p">,</span>
                 <span class="n">initial_values</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Returns a vector-autoregressive process with correlated innovations.</span>

<span class="sd">    Useful for testing.</span>

<span class="sd">    Example:</span>
<span class="sd">        graph=numpy.array([[[0.2,0.,0.],[0.5,0.,0.]],</span>
<span class="sd">                           [[0.,0.1,0. ],[0.3,0.,0.]]])</span>

<span class="sd">        represents a process</span>

<span class="sd">        X_1(t) = 0.2 X_1(t-1) + 0.5 X_2(t-1) + eps_1(t)</span>
<span class="sd">        X_2(t) = 0.3 X_2(t-1) + 0.1 X_1(t-2) + eps_2(t)</span>

<span class="sd">        with inv_inno_cov being the negative (except for diagonal) inverse</span>
<span class="sd">        covariance matrix of (eps_1(t), eps_2(t)) OR inno_cov being</span>
<span class="sd">        the covariance. Initial values can also be provided.</span>


<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    graph : array</span>
<span class="sd">        Lagged connectivity matrices. Shape is (n_nodes, n_nodes, max_delay+1)</span>
<span class="sd">    add_noise : bool, optional (default: True)</span>
<span class="sd">        Flag to add random noise or not</span>
<span class="sd">    inno_cov : array, optional (default: None)</span>
<span class="sd">        Covariance matrix of innovations.</span>
<span class="sd">    invert_inno : bool, optional (defualt : False)</span>
<span class="sd">        Flag to negate off-diagonal elements of inno_cov and invert it before</span>
<span class="sd">        using it as the covariance matrix of innovations</span>
<span class="sd">    T : int, optional (default: 100)</span>
<span class="sd">        Sample size.</span>

<span class="sd">    initial_values : array, optional (defult: None)</span>
<span class="sd">        Initial values for each node. Shape is (n_nodes, max_delay+1), i.e. must</span>
<span class="sd">        be of shape (graph.shape[1], graph.shape[2]).</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    X : array</span>
<span class="sd">        Array of realization.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">n_nodes</span><span class="p">,</span> <span class="n">_</span><span class="p">,</span> <span class="n">period</span> <span class="o">=</span> <span class="n">graph</span><span class="o">.</span><span class="n">shape</span>

    <span class="n">time</span> <span class="o">=</span> <span class="n">T</span>
    <span class="c1"># Test stability</span>
    <span class="n">_check_stability</span><span class="p">(</span><span class="n">graph</span><span class="p">)</span>

    <span class="c1"># Generate the returned data</span>
    <span class="n">data</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">randn</span><span class="p">(</span><span class="n">n_nodes</span><span class="p">,</span> <span class="n">time</span><span class="p">)</span>
    <span class="c1"># Load the initial values</span>
    <span class="k">if</span> <span class="n">initial_values</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
        <span class="c1"># Check the shape of the initial values</span>
        <span class="n">_check_initial_values</span><span class="p">(</span><span class="n">initial_values</span><span class="p">,</span> <span class="n">data</span><span class="p">[:,</span> <span class="p">:</span><span class="n">period</span><span class="p">]</span><span class="o">.</span><span class="n">shape</span><span class="p">)</span>
        <span class="c1"># Input the initial values</span>
        <span class="n">data</span><span class="p">[:,</span> <span class="p">:</span><span class="n">period</span><span class="p">]</span> <span class="o">=</span> <span class="n">initial_values</span>

    <span class="c1"># Check if we are adding noise</span>
    <span class="n">noise</span> <span class="o">=</span> <span class="kc">None</span>
    <span class="k">if</span> <span class="n">add_noise</span><span class="p">:</span>
        <span class="c1"># Use inno_cov if it was provided</span>
        <span class="k">if</span> <span class="n">inno_cov</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
            <span class="n">noise</span> <span class="o">=</span> <span class="n">_generate_noise</span><span class="p">(</span><span class="n">inno_cov</span><span class="p">,</span>
                                    <span class="n">time</span><span class="o">=</span><span class="n">time</span><span class="p">,</span>
                                    <span class="n">use_inverse</span><span class="o">=</span><span class="n">invert_inno</span><span class="p">)</span>
        <span class="c1"># Otherwise just use uncorrelated random noise</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">noise</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">randn</span><span class="p">(</span><span class="n">time</span><span class="p">,</span> <span class="n">n_nodes</span><span class="p">)</span>

    <span class="k">for</span> <span class="n">a_time</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">period</span><span class="p">,</span> <span class="n">time</span><span class="p">):</span>
        <span class="n">data_past</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">repeat</span><span class="p">(</span>
            <span class="n">data</span><span class="p">[:,</span> <span class="n">a_time</span><span class="o">-</span><span class="n">period</span><span class="p">:</span><span class="n">a_time</span><span class="p">][:,</span> <span class="p">::</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="n">n_nodes</span><span class="p">,</span> <span class="n">period</span><span class="p">),</span>
            <span class="n">n_nodes</span><span class="p">,</span> <span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
        <span class="n">data</span><span class="p">[:,</span> <span class="n">a_time</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">data_past</span><span class="o">*</span><span class="n">graph</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">2</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
        <span class="k">if</span> <span class="n">add_noise</span><span class="p">:</span>
            <span class="n">data</span><span class="p">[:,</span> <span class="n">a_time</span><span class="p">]</span> <span class="o">+=</span> <span class="n">noise</span><span class="p">[</span><span class="n">a_time</span><span class="p">]</span>

    <span class="k">return</span> <span class="n">data</span><span class="o">.</span><span class="n">transpose</span><span class="p">()</span>

<span class="k">def</span> <span class="nf">_iter_coeffs</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Iterator through the current parents_neighbors_coeffs structure.  Mainly to</span>
<span class="sd">    save repeated code and make it easier to change this structure.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    parents_neighbors_coeffs : dict</span>
<span class="sd">        Dictionary of format:</span>
<span class="sd">        {..., j:[((var1, lag1), coef1), ((var2, lag2), coef2), ...], ...} for</span>
<span class="sd">        all variables where vars must be in [0..N-1] and lags &lt;= 0 with number</span>
<span class="sd">        of variables N.</span>

<span class="sd">    Yields</span>
<span class="sd">    -------</span>
<span class="sd">    (node_id, parent_id, time_lag, coeff) : tuple</span>
<span class="sd">        Tuple defining the relationship between nodes across time</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Iterate through all defined nodes</span>
    <span class="k">for</span> <span class="n">node_id</span> <span class="ow">in</span> <span class="nb">list</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
        <span class="c1"># Iterate over parent nodes and unpack node and coeff</span>
        <span class="k">for</span> <span class="p">(</span><span class="n">parent_id</span><span class="p">,</span> <span class="n">time_lag</span><span class="p">),</span> <span class="n">coeff</span> <span class="ow">in</span> <span class="n">parents_neighbors_coeffs</span><span class="p">[</span><span class="n">node_id</span><span class="p">]:</span>
            <span class="c1"># Yield the entry</span>
            <span class="k">yield</span> <span class="n">node_id</span><span class="p">,</span> <span class="n">parent_id</span><span class="p">,</span> <span class="n">time_lag</span><span class="p">,</span> <span class="n">coeff</span>

<span class="k">def</span> <span class="nf">_check_parent_neighbor</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Checks to insure input parent-neighbor connectivity input is sane.  This</span>
<span class="sd">    means that:</span>
<span class="sd">        * all time lags are non-positive</span>
<span class="sd">        * all parent nodes are included as nodes themselves</span>
<span class="sd">        * all node indexing is contiguous</span>
<span class="sd">        * all node indexing starts from zero</span>
<span class="sd">    Raises a ValueError if any one of these conditions are not met.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    parents_neighbors_coeffs : dict</span>
<span class="sd">        Dictionary of format:</span>
<span class="sd">        {..., j:[((var1, lag1), coef1), ((var2, lag2), coef2), ...], ...} for</span>
<span class="sd">        all variables where vars must be in [0..N-1] and lags &lt;= 0 with number</span>
<span class="sd">        of variables N.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Initialize some lists for checking later</span>
    <span class="n">all_nodes</span> <span class="o">=</span> <span class="nb">set</span><span class="p">()</span>
    <span class="n">all_parents</span> <span class="o">=</span> <span class="nb">set</span><span class="p">()</span>
    <span class="c1"># Iterate through variables</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">list</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
        <span class="c1"># Cache all node ids to ensure they are contiguous</span>
        <span class="n">all_nodes</span><span class="o">.</span><span class="n">add</span><span class="p">(</span><span class="n">j</span><span class="p">)</span>
    <span class="c1"># Iterate through all nodes</span>
    <span class="k">for</span> <span class="n">j</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">tau</span><span class="p">,</span> <span class="n">_</span> <span class="ow">in</span> <span class="n">_iter_coeffs</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
        <span class="c1"># Check all time lags are equal to or less than zero</span>
        <span class="k">if</span> <span class="n">tau</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;Lag between parent </span><span class="si">{}</span><span class="s2"> and node </span><span class="si">{}</span><span class="s2">&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">j</span><span class="p">)</span><span class="o">+</span>\
                             <span class="s2">&quot; is </span><span class="si">{}</span><span class="s2"> &gt; 0, must be &lt;= 0!&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">tau</span><span class="p">))</span>
        <span class="c1"># Cache all parent ids to ensure they are mentioned as node ids</span>
        <span class="n">all_parents</span><span class="o">.</span><span class="n">add</span><span class="p">(</span><span class="n">i</span><span class="p">)</span>
    <span class="c1"># Check that all nodes are contiguous from zero</span>
    <span class="n">all_nodes_list</span> <span class="o">=</span> <span class="nb">sorted</span><span class="p">(</span><span class="nb">list</span><span class="p">(</span><span class="n">all_nodes</span><span class="p">))</span>
    <span class="k">if</span> <span class="n">all_nodes_list</span> <span class="o">!=</span> <span class="nb">list</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">all_nodes_list</span><span class="p">))):</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;Node IDs in input dictionary must be contiguous&quot;</span><span class="o">+</span>\
                         <span class="s2">&quot; and start from zero!</span><span class="se">\n</span><span class="s2">&quot;</span><span class="o">+</span>\
                         <span class="s2">&quot; Found IDs : [&quot;</span> <span class="o">+</span>\
                         <span class="s2">&quot;,&quot;</span><span class="o">.</span><span class="n">join</span><span class="p">(</span><span class="nb">map</span><span class="p">(</span><span class="nb">str</span><span class="p">,</span> <span class="n">all_nodes_list</span><span class="p">))</span><span class="o">+</span> <span class="s2">&quot;]&quot;</span><span class="p">)</span>
    <span class="c1"># Check that all parent nodes are mentioned as a node ID</span>
    <span class="k">if</span> <span class="ow">not</span> <span class="n">all_parents</span><span class="o">.</span><span class="n">issubset</span><span class="p">(</span><span class="n">all_nodes</span><span class="p">):</span>
        <span class="n">missing_nodes</span> <span class="o">=</span> <span class="nb">sorted</span><span class="p">(</span><span class="nb">list</span><span class="p">(</span><span class="n">all_parents</span> <span class="o">-</span> <span class="n">all_nodes</span><span class="p">))</span>
        <span class="n">all_parents_list</span> <span class="o">=</span> <span class="nb">sorted</span><span class="p">(</span><span class="nb">list</span><span class="p">(</span><span class="n">all_parents</span><span class="p">))</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;Parent IDs in input dictionary must also be in set&quot;</span><span class="o">+</span>\
                         <span class="s2">&quot; of node IDs.&quot;</span><span class="o">+</span>\
                         <span class="s2">&quot;</span><span class="se">\n</span><span class="s2"> Parent IDs &quot;</span><span class="o">+</span><span class="s2">&quot; &quot;</span><span class="o">.</span><span class="n">join</span><span class="p">(</span><span class="nb">map</span><span class="p">(</span><span class="nb">str</span><span class="p">,</span> <span class="n">all_parents_list</span><span class="p">))</span><span class="o">+</span>\
                         <span class="s2">&quot;</span><span class="se">\n</span><span class="s2"> Node IDs &quot;</span><span class="o">+</span><span class="s2">&quot; &quot;</span><span class="o">.</span><span class="n">join</span><span class="p">(</span><span class="nb">map</span><span class="p">(</span><span class="nb">str</span><span class="p">,</span> <span class="n">all_nodes_list</span><span class="p">))</span> <span class="o">+</span>\
                         <span class="s2">&quot;</span><span class="se">\n</span><span class="s2"> Missing IDs &quot;</span> <span class="o">+</span> <span class="s2">&quot; &quot;</span><span class="o">.</span><span class="n">join</span><span class="p">(</span><span class="nb">map</span><span class="p">(</span><span class="nb">str</span><span class="p">,</span> <span class="n">missing_nodes</span><span class="p">)))</span>

<span class="k">def</span> <span class="nf">_check_symmetric_relations</span><span class="p">(</span><span class="n">a_matrix</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Check if the argument matrix is symmetric.  Raise a value error with details</span>
<span class="sd">    about the offending elements if it is not.  This is useful for checking the</span>
<span class="sd">    instantaneously linked nodes have the same link strength.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    a_matrix : 2D numpy array</span>
<span class="sd">        Relationships between nodes at tau = 0. Indexed such that first index is</span>
<span class="sd">        node and second is parent, i.e. node j with parent i has strength</span>
<span class="sd">        a_matrix[j,i]</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Check it is symmetric</span>
    <span class="k">if</span> <span class="ow">not</span> <span class="n">np</span><span class="o">.</span><span class="n">allclose</span><span class="p">(</span><span class="n">a_matrix</span><span class="p">,</span> <span class="n">a_matrix</span><span class="o">.</span><span class="n">T</span><span class="p">,</span> <span class="n">rtol</span><span class="o">=</span><span class="mf">1e-10</span><span class="p">,</span> <span class="n">atol</span><span class="o">=</span><span class="mf">1e-10</span><span class="p">):</span>
        <span class="c1"># Store the disagreement elements</span>
        <span class="n">bad_elems</span> <span class="o">=</span> <span class="o">~</span><span class="n">np</span><span class="o">.</span><span class="n">isclose</span><span class="p">(</span><span class="n">a_matrix</span><span class="p">,</span> <span class="n">a_matrix</span><span class="o">.</span><span class="n">T</span><span class="p">,</span> <span class="n">rtol</span><span class="o">=</span><span class="mf">1e-10</span><span class="p">,</span> <span class="n">atol</span><span class="o">=</span><span class="mf">1e-10</span><span class="p">)</span>
        <span class="n">bad_idxs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">argwhere</span><span class="p">(</span><span class="n">bad_elems</span><span class="p">)</span>
        <span class="n">error_message</span> <span class="o">=</span> <span class="s2">&quot;&quot;</span>
        <span class="k">for</span> <span class="n">node</span><span class="p">,</span> <span class="n">parent</span> <span class="ow">in</span> <span class="n">bad_idxs</span><span class="p">:</span>
            <span class="c1"># Check that we haven&#39;t already printed about this pair</span>
            <span class="k">if</span> <span class="n">bad_elems</span><span class="p">[</span><span class="n">node</span><span class="p">,</span> <span class="n">parent</span><span class="p">]:</span>
                <span class="n">error_message</span> <span class="o">+=</span> \
                    <span class="s2">&quot;Parent </span><span class="si">{:d}</span><span class="s2"> of node </span><span class="si">{:d}</span><span class="s2">&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">parent</span><span class="p">,</span> <span class="n">node</span><span class="p">)</span><span class="o">+</span>\
                    <span class="s2">&quot; has coefficient </span><span class="si">{:f}</span><span class="s2">.</span><span class="se">\n</span><span class="s2">&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">a_matrix</span><span class="p">[</span><span class="n">node</span><span class="p">,</span> <span class="n">parent</span><span class="p">])</span><span class="o">+</span>\
                    <span class="s2">&quot;Parent </span><span class="si">{:d}</span><span class="s2"> of node </span><span class="si">{:d}</span><span class="s2">&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">node</span><span class="p">,</span> <span class="n">parent</span><span class="p">)</span><span class="o">+</span>\
                    <span class="s2">&quot; has coefficient </span><span class="si">{:f}</span><span class="s2">.</span><span class="se">\n</span><span class="s2">&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">a_matrix</span><span class="p">[</span><span class="n">parent</span><span class="p">,</span> <span class="n">node</span><span class="p">])</span>
            <span class="c1"># Check if we already printed about this one</span>
            <span class="n">bad_elems</span><span class="p">[</span><span class="n">node</span><span class="p">,</span> <span class="n">parent</span><span class="p">]</span> <span class="o">=</span> <span class="kc">False</span>
            <span class="n">bad_elems</span><span class="p">[</span><span class="n">parent</span><span class="p">,</span> <span class="n">node</span><span class="p">]</span> <span class="o">=</span> <span class="kc">False</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;Relationships between nodes at tau=0 are not&quot;</span><span class="o">+</span>\
                         <span class="s2">&quot; symmetric!</span><span class="se">\n</span><span class="s2">&quot;</span><span class="o">+</span><span class="n">error_message</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">_find_max_time_lag_and_node_id</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Function to find the maximum time lag in the parent-neighbors-coefficients</span>
<span class="sd">    object, as well as the largest node ID</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    parents_neighbors_coeffs : dict</span>
<span class="sd">        Dictionary of format:</span>
<span class="sd">        {..., j:[((var1, lag1), coef1), ((var2, lag2), coef2), ...], ...} for</span>
<span class="sd">        all variables where vars must be in [0..N-1] and lags &lt;= 0 with number</span>
<span class="sd">        of variables N.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    (max_time_lag, max_node_id) : tuple</span>
<span class="sd">        Tuple of the maximum time lag and maximum node ID</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Default maximum lag and node ID</span>
    <span class="n">max_time_lag</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">max_node_id</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="o">.</span><span class="n">keys</span><span class="p">())</span> <span class="o">-</span> <span class="mi">1</span>
    <span class="c1"># Iterate through the keys in parents_neighbors_coeffs</span>
    <span class="k">for</span> <span class="n">j</span><span class="p">,</span> <span class="n">_</span><span class="p">,</span> <span class="n">tau</span><span class="p">,</span> <span class="n">_</span> <span class="ow">in</span> <span class="n">_iter_coeffs</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
        <span class="c1"># Find max lag time</span>
        <span class="n">max_time_lag</span> <span class="o">=</span> <span class="nb">max</span><span class="p">(</span><span class="n">max_time_lag</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">tau</span><span class="p">))</span>
        <span class="c1"># Find the max node ID</span>
        <span class="c1"># max_node_id = max(max_node_id, j)</span>
    <span class="c1"># Return these values</span>
    <span class="k">return</span> <span class="n">max_time_lag</span><span class="p">,</span> <span class="n">max_node_id</span>

<span class="k">def</span> <span class="nf">_get_true_parent_neighbor_dict</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Function to return the dictionary of true parent neighbor causal</span>
<span class="sd">    connections in time.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    parents_neighbors_coeffs : dict</span>
<span class="sd">        Dictionary of format:</span>
<span class="sd">        {..., j:[((var1, lag1), coef1), ((var2, lag2), coef2), ...], ...} for</span>
<span class="sd">        all variables where vars must be in [0..N-1] and lags &lt;= 0 with number</span>
<span class="sd">        of variables N.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    true_parent_neighbor : dict</span>
<span class="sd">        Dictionary of lists of tuples.  The dictionary is keyed by node ID, the</span>
<span class="sd">        list stores the tuple values (parent_node_id, time_lag)</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Initialize the returned dictionary of lists</span>
    <span class="n">true_parents_neighbors</span> <span class="o">=</span> <span class="n">defaultdict</span><span class="p">(</span><span class="nb">list</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">parents_neighbors_coeffs</span><span class="p">:</span>
        <span class="k">for</span> <span class="n">link_props</span> <span class="ow">in</span> <span class="n">parents_neighbors_coeffs</span><span class="p">[</span><span class="n">j</span><span class="p">]:</span>
            <span class="n">i</span><span class="p">,</span> <span class="n">tau</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
            <span class="n">coeff</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
            <span class="c1"># Add parent node id and lag if non-zero coeff</span>
            <span class="k">if</span> <span class="n">coeff</span> <span class="o">!=</span> <span class="mf">0.</span><span class="p">:</span>
                <span class="n">true_parents_neighbors</span><span class="p">[</span><span class="n">j</span><span class="p">]</span><span class="o">.</span><span class="n">append</span><span class="p">((</span><span class="n">i</span><span class="p">,</span> <span class="n">tau</span><span class="p">))</span>
    <span class="c1"># Return the true relations</span>
    <span class="k">return</span> <span class="n">true_parents_neighbors</span>

<span class="k">def</span> <span class="nf">_get_covariance_matrix</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Determines the covariance matrix for correlated innovations</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    parents_neighbors_coeffs : dict</span>
<span class="sd">        Dictionary of format:</span>
<span class="sd">        {..., j:[((var1, lag1), coef1), ((var2, lag2), coef2), ...], ...} for</span>
<span class="sd">        all variables where vars must be in [0..N-1] and lags &lt;= 0 with number</span>
<span class="sd">        of variables N.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    covar_matrix : numpy array</span>
<span class="sd">        Covariance matrix implied by the parents_neighbors_coeffs.  Used to</span>
<span class="sd">        generate correlated innovations.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the total number of nodes</span>
    <span class="n">_</span><span class="p">,</span> <span class="n">max_node_id</span> <span class="o">=</span> \
            <span class="n">_find_max_time_lag_and_node_id</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">)</span>
    <span class="n">n_nodes</span> <span class="o">=</span> <span class="n">max_node_id</span> <span class="o">+</span> <span class="mi">1</span>
    <span class="c1"># Initialize the covariance matrix</span>
    <span class="n">covar_matrix</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">identity</span><span class="p">(</span><span class="n">n_nodes</span><span class="p">)</span>
    <span class="c1"># Iterate through all the node connections</span>
    <span class="k">for</span> <span class="n">j</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">tau</span><span class="p">,</span> <span class="n">coeff</span> <span class="ow">in</span> <span class="n">_iter_coeffs</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
        <span class="c1"># Add to covar_matrix if node connection is instantaneous</span>
        <span class="k">if</span> <span class="n">tau</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">covar_matrix</span><span class="p">[</span><span class="n">j</span><span class="p">,</span> <span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">coeff</span>
    <span class="k">return</span> <span class="n">covar_matrix</span>

<span class="k">def</span> <span class="nf">_get_lag_connect_matrix</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Generates the lagged connectivity matrix from a parent-neighbor</span>
<span class="sd">    connectivity dictionary.  Used to generate the input for _var_network</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    parents_neighbors_coeffs : dict</span>
<span class="sd">        Dictionary of format:</span>
<span class="sd">        {..., j:[((var1, lag1), coef1), ((var2, lag2), coef2), ...], ...} for</span>
<span class="sd">        all variables where vars must be in [0..N-1] and lags &lt;= 0 with number</span>
<span class="sd">        of variables N.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    connect_matrix : numpy array</span>
<span class="sd">        Lagged connectivity matrix. Shape is (n_nodes, n_nodes, max_delay+1)</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the total number of nodes and time lag</span>
    <span class="n">max_time_lag</span><span class="p">,</span> <span class="n">max_node_id</span> <span class="o">=</span> \
            <span class="n">_find_max_time_lag_and_node_id</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">)</span>
    <span class="n">n_nodes</span> <span class="o">=</span> <span class="n">max_node_id</span> <span class="o">+</span> <span class="mi">1</span>
    <span class="n">n_times</span> <span class="o">=</span> <span class="n">max_time_lag</span> <span class="o">+</span> <span class="mi">1</span>
    <span class="c1"># Initialize full time graph</span>
    <span class="n">connect_matrix</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">n_nodes</span><span class="p">,</span> <span class="n">n_nodes</span><span class="p">,</span> <span class="n">n_times</span><span class="p">))</span>
    <span class="k">for</span> <span class="n">j</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">tau</span><span class="p">,</span> <span class="n">coeff</span> <span class="ow">in</span> <span class="n">_iter_coeffs</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">):</span>
        <span class="c1"># If there is a non-zero time lag, add the connection to the matrix</span>
        <span class="k">if</span> <span class="n">tau</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">connect_matrix</span><span class="p">[</span><span class="n">j</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="o">-</span><span class="p">(</span><span class="n">tau</span><span class="o">+</span><span class="mi">1</span><span class="p">)]</span> <span class="o">=</span> <span class="n">coeff</span>
    <span class="c1"># Return the connectivity matrix</span>
    <span class="k">return</span> <span class="n">connect_matrix</span>

<div class="viewcode-block" id="var_process"><a class="viewcode-back" href="../../../index.html#tigramite.toymodels.structural_causal_processes.var_process">[docs]</a><span class="k">def</span> <span class="nf">var_process</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">,</span> <span class="n">T</span><span class="o">=</span><span class="mi">1000</span><span class="p">,</span> <span class="n">use</span><span class="o">=</span><span class="s1">&#39;inv_inno_cov&#39;</span><span class="p">,</span>
                <span class="n">verbosity</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">initial_values</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;Returns a vector-autoregressive process with correlated innovations.</span>

<span class="sd">    Wrapper around var_network with possibly more user-friendly input options.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    parents_neighbors_coeffs : dict</span>
<span class="sd">        Dictionary of format: {..., j:[((var1, lag1), coef1), ((var2, lag2),</span>
<span class="sd">        coef2), ...], ...} for all variables where vars must be in [0..N-1]</span>
<span class="sd">        and lags &lt;= 0 with number of variables N. If lag=0, a nonzero value</span>
<span class="sd">        in the covariance matrix (or its inverse) is implied. These should be</span>
<span class="sd">        the same for (i, j) and (j, i).</span>
<span class="sd">    use : str, optional (default: &#39;inv_inno_cov&#39;)</span>
<span class="sd">        Specifier, either &#39;inno_cov&#39; or &#39;inv_inno_cov&#39;.</span>
<span class="sd">        Any other specifier will result in non-correlated noise.</span>
<span class="sd">        For debugging, &#39;no_noise&#39; can also be specified, in which case random</span>
<span class="sd">        noise will be disabled.</span>
<span class="sd">    T : int, optional (default: 1000)</span>
<span class="sd">        Sample size.</span>
<span class="sd">    verbosity : int, optional (default: 0)</span>
<span class="sd">        Level of verbosity.</span>
<span class="sd">    initial_values : array, optional (default: None)</span>
<span class="sd">        Initial values for each node. Shape must be (N, max_delay+1)</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    data : array-like</span>
<span class="sd">        Data generated from this process</span>
<span class="sd">    true_parent_neighbor : dict</span>
<span class="sd">        Dictionary of lists of tuples.  The dictionary is keyed by node ID, the</span>
<span class="sd">        list stores the tuple values (parent_node_id, time_lag)</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Check the input parents_neighbors_coeffs dictionary for sanity</span>
    <span class="n">_check_parent_neighbor</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">)</span>
    <span class="c1"># Generate the true parent neighbors graph</span>
    <span class="n">true_parents_neighbors</span> <span class="o">=</span> \
        <span class="n">_get_true_parent_neighbor_dict</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">)</span>
    <span class="c1"># Generate the correlated innovations</span>
    <span class="n">innos</span> <span class="o">=</span> <span class="n">_get_covariance_matrix</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">)</span>
    <span class="c1"># Generate the lagged connectivity matrix for _var_network</span>
    <span class="n">connect_matrix</span> <span class="o">=</span> <span class="n">_get_lag_connect_matrix</span><span class="p">(</span><span class="n">parents_neighbors_coeffs</span><span class="p">)</span>
    <span class="c1"># Default values as per &#39;inno_cov&#39;</span>
    <span class="n">add_noise</span> <span class="o">=</span> <span class="kc">True</span>
    <span class="n">invert_inno</span> <span class="o">=</span> <span class="kc">False</span>
    <span class="c1"># Use the correlated innovations</span>
    <span class="k">if</span> <span class="n">use</span> <span class="o">==</span> <span class="s1">&#39;inno_cov&#39;</span><span class="p">:</span>
        <span class="k">if</span> <span class="n">verbosity</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="nb">print</span><span class="p">(</span><span class="s2">&quot;</span><span class="se">\n</span><span class="s2">Innovation Cov =</span><span class="se">\n</span><span class="si">%s</span><span class="s2">&quot;</span> <span class="o">%</span> <span class="nb">str</span><span class="p">(</span><span class="n">innos</span><span class="p">))</span>
    <span class="c1"># Use the inverted correlated innovations</span>
    <span class="k">elif</span> <span class="n">use</span> <span class="o">==</span> <span class="s1">&#39;inv_inno_cov&#39;</span><span class="p">:</span>
        <span class="n">invert_inno</span> <span class="o">=</span> <span class="kc">True</span>
        <span class="k">if</span> <span class="n">verbosity</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="nb">print</span><span class="p">(</span><span class="s2">&quot;</span><span class="se">\n</span><span class="s2">Inverse Innovation Cov =</span><span class="se">\n</span><span class="si">%s</span><span class="s2">&quot;</span> <span class="o">%</span> <span class="nb">str</span><span class="p">(</span><span class="n">innos</span><span class="p">))</span>
    <span class="c1"># Do not use any noise</span>
    <span class="k">elif</span> <span class="n">use</span> <span class="o">==</span> <span class="s1">&#39;no_noise&#39;</span><span class="p">:</span>
        <span class="n">add_noise</span> <span class="o">=</span> <span class="kc">False</span>
        <span class="k">if</span> <span class="n">verbosity</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="nb">print</span><span class="p">(</span><span class="s2">&quot;</span><span class="se">\n</span><span class="s2">Inverse Innovation Cov =</span><span class="se">\n</span><span class="si">%s</span><span class="s2">&quot;</span> <span class="o">%</span> <span class="nb">str</span><span class="p">(</span><span class="n">innos</span><span class="p">))</span>
    <span class="c1"># Use decorrelated noise</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">innos</span> <span class="o">=</span> <span class="kc">None</span>
    <span class="c1"># Ensure the innovation matrix is symmetric if it is used</span>
    <span class="k">if</span> <span class="p">(</span><span class="n">innos</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">)</span> <span class="ow">and</span> <span class="n">add_noise</span><span class="p">:</span>
        <span class="n">_check_symmetric_relations</span><span class="p">(</span><span class="n">innos</span><span class="p">)</span>
    <span class="c1"># Generate the data using _var_network</span>
    <span class="n">data</span> <span class="o">=</span> <span class="n">_var_network</span><span class="p">(</span><span class="n">graph</span><span class="o">=</span><span class="n">connect_matrix</span><span class="p">,</span>
                        <span class="n">add_noise</span><span class="o">=</span><span class="n">add_noise</span><span class="p">,</span>
                        <span class="n">inno_cov</span><span class="o">=</span><span class="n">innos</span><span class="p">,</span>
                        <span class="n">invert_inno</span><span class="o">=</span><span class="n">invert_inno</span><span class="p">,</span>
                        <span class="n">T</span><span class="o">=</span><span class="n">T</span><span class="p">,</span>
                        <span class="n">initial_values</span><span class="o">=</span><span class="n">initial_values</span><span class="p">)</span>
    <span class="c1"># Return the data</span>
    <span class="k">return</span> <span class="n">data</span><span class="p">,</span> <span class="n">true_parents_neighbors</span></div>

<span class="k">class</span> <span class="nc">_Graph</span><span class="p">():</span>
<span class="w">    </span><span class="sa">r</span><span class="sd">&quot;&quot;&quot;Helper class to handle graph properties.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    vertices : list</span>
<span class="sd">        List of nodes.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">def</span> <span class="fm">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span><span class="n">vertices</span><span class="p">):</span> 
        <span class="bp">self</span><span class="o">.</span><span class="n">graph</span> <span class="o">=</span> <span class="n">defaultdict</span><span class="p">(</span><span class="nb">list</span><span class="p">)</span> 
        <span class="bp">self</span><span class="o">.</span><span class="n">V</span> <span class="o">=</span> <span class="n">vertices</span> 
  
    <span class="k">def</span> <span class="nf">addEdge</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span><span class="n">u</span><span class="p">,</span><span class="n">v</span><span class="p">):</span>
<span class="w">        </span><span class="sd">&quot;&quot;&quot;Adding edge to graph.&quot;&quot;&quot;</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">graph</span><span class="p">[</span><span class="n">u</span><span class="p">]</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">v</span><span class="p">)</span> 
  
    <span class="k">def</span> <span class="nf">isCyclicUtil</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">v</span><span class="p">,</span> <span class="n">visited</span><span class="p">,</span> <span class="n">recStack</span><span class="p">):</span> 
<span class="w">        </span><span class="sd">&quot;&quot;&quot;Utility function to return whether graph is cyclic.&quot;&quot;&quot;</span>
        <span class="c1"># Mark current node as visited and</span>
        <span class="c1"># adds to recursion stack </span>
        <span class="n">visited</span><span class="p">[</span><span class="n">v</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>
        <span class="n">recStack</span><span class="p">[</span><span class="n">v</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>
  
        <span class="c1"># Recur for all neighbours </span>
        <span class="c1"># if any neighbour is visited and in  </span>
        <span class="c1"># recStack then graph is cyclic </span>
        <span class="k">for</span> <span class="n">neighbour</span> <span class="ow">in</span> <span class="bp">self</span><span class="o">.</span><span class="n">graph</span><span class="p">[</span><span class="n">v</span><span class="p">]:</span> 
            <span class="k">if</span> <span class="n">visited</span><span class="p">[</span><span class="n">neighbour</span><span class="p">]</span> <span class="o">==</span> <span class="kc">False</span><span class="p">:</span> 
                <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">isCyclicUtil</span><span class="p">(</span><span class="n">neighbour</span><span class="p">,</span> <span class="n">visited</span><span class="p">,</span> <span class="n">recStack</span><span class="p">)</span> <span class="o">==</span> <span class="kc">True</span><span class="p">:</span> 
                    <span class="k">return</span> <span class="kc">True</span>
            <span class="k">elif</span> <span class="n">recStack</span><span class="p">[</span><span class="n">neighbour</span><span class="p">]</span> <span class="o">==</span> <span class="kc">True</span><span class="p">:</span> 
                <span class="k">return</span> <span class="kc">True</span>
  
        <span class="c1"># The node needs to be poped from  </span>
        <span class="c1"># recursion stack before function ends </span>
        <span class="n">recStack</span><span class="p">[</span><span class="n">v</span><span class="p">]</span> <span class="o">=</span> <span class="kc">False</span>
        <span class="k">return</span> <span class="kc">False</span>
  
    <span class="k">def</span> <span class="nf">isCyclic</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
<span class="w">        </span><span class="sd">&quot;&quot;&quot;Returns whether graph is cyclic.&quot;&quot;&quot;</span>
        <span class="n">visited</span> <span class="o">=</span> <span class="p">[</span><span class="kc">False</span><span class="p">]</span> <span class="o">*</span> <span class="bp">self</span><span class="o">.</span><span class="n">V</span> 
        <span class="n">recStack</span> <span class="o">=</span> <span class="p">[</span><span class="kc">False</span><span class="p">]</span> <span class="o">*</span> <span class="bp">self</span><span class="o">.</span><span class="n">V</span> 
        <span class="k">for</span> <span class="n">node</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">V</span><span class="p">):</span> 
            <span class="k">if</span> <span class="n">visited</span><span class="p">[</span><span class="n">node</span><span class="p">]</span> <span class="o">==</span> <span class="kc">False</span><span class="p">:</span> 
                <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">isCyclicUtil</span><span class="p">(</span><span class="n">node</span><span class="p">,</span><span class="n">visited</span><span class="p">,</span><span class="n">recStack</span><span class="p">)</span> <span class="o">==</span> <span class="kc">True</span><span class="p">:</span> 
                    <span class="k">return</span> <span class="kc">True</span>
        <span class="k">return</span> <span class="kc">False</span>
  
    <span class="k">def</span> <span class="nf">topologicalSortUtil</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span><span class="n">v</span><span class="p">,</span><span class="n">visited</span><span class="p">,</span><span class="n">stack</span><span class="p">):</span>
<span class="w">        </span><span class="sd">&quot;&quot;&quot;A recursive function used by topologicalSort .&quot;&quot;&quot;</span>
        <span class="c1"># Mark the current node as visited.</span>
        <span class="n">visited</span><span class="p">[</span><span class="n">v</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>

        <span class="c1"># Recur for all the vertices adjacent to this vertex</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="bp">self</span><span class="o">.</span><span class="n">graph</span><span class="p">[</span><span class="n">v</span><span class="p">]:</span>
            <span class="k">if</span> <span class="n">visited</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">==</span> <span class="kc">False</span><span class="p">:</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">topologicalSortUtil</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">visited</span><span class="p">,</span><span class="n">stack</span><span class="p">)</span>

        <span class="c1"># Push current vertex to stack which stores result</span>
        <span class="n">stack</span><span class="o">.</span><span class="n">insert</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="n">v</span><span class="p">)</span>

    <span class="k">def</span> <span class="nf">topologicalSort</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
<span class="w">        </span><span class="sd">&quot;&quot;&quot;A sorting function. &quot;&quot;&quot;</span>
        <span class="c1"># Mark all the vertices as not visited </span>
        <span class="n">visited</span> <span class="o">=</span> <span class="p">[</span><span class="kc">False</span><span class="p">]</span><span class="o">*</span><span class="bp">self</span><span class="o">.</span><span class="n">V</span> 
        <span class="n">stack</span> <span class="o">=</span><span class="p">[]</span> 

        <span class="c1"># Call the recursive helper function to store Topological </span>
        <span class="c1"># Sort starting from all vertices one by one </span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">V</span><span class="p">):</span> 
          <span class="k">if</span> <span class="n">visited</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">==</span> <span class="kc">False</span><span class="p">:</span> 
              <span class="bp">self</span><span class="o">.</span><span class="n">topologicalSortUtil</span><span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">visited</span><span class="p">,</span><span class="n">stack</span><span class="p">)</span> 

        <span class="k">return</span> <span class="n">stack</span>

<div class="viewcode-block" id="structural_causal_process"><a class="viewcode-back" href="../../../index.html#tigramite.toymodels.structural_causal_processes.structural_causal_process">[docs]</a><span class="k">def</span> <span class="nf">structural_causal_process</span><span class="p">(</span><span class="n">links</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="n">noises</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> 
                        <span class="n">intervention</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">intervention_type</span><span class="o">=</span><span class="s1">&#39;hard&#39;</span><span class="p">,</span>
                        <span class="n">transient_fraction</span><span class="o">=</span><span class="mf">0.2</span><span class="p">,</span>
                        <span class="n">seed</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;Returns a time series generated from a structural causal process.</span>

<span class="sd">    Allows lagged and contemporaneous dependencies and includes the option</span>
<span class="sd">    to have intervened variables or particular samples.</span>

<span class="sd">    The interventional data is in particular useful for generating ground</span>
<span class="sd">    truth for the CausalEffects class.</span>

<span class="sd">    In more detail, the method implements a generalized additive noise model process of the form</span>

<span class="sd">    .. math:: X^j_t = \\eta^j_t + \\sum_{X^i_{t-\\tau}\\in \\mathcal{P}(X^j_t)}</span>
<span class="sd">              c^i_{\\tau} f^i_{\\tau}(X^i_{t-\\tau})</span>

<span class="sd">    Links have the format ``{0:[((i, -tau), coeff, func),...], 1:[...],</span>
<span class="sd">    ...}`` where ``func`` can be an arbitrary (nonlinear) function provided</span>
<span class="sd">    as a python callable with one argument and coeff is the multiplication</span>
<span class="sd">    factor. The noise distributions of :math:`\\eta^j` can be specified in</span>
<span class="sd">    ``noises``.</span>

<span class="sd">    Through the parameters ``intervention`` and ``intervention_type`` the model</span>
<span class="sd">    can also be generated with intervened variables.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    links : dict</span>
<span class="sd">        Dictionary of format: {0:[((i, -tau), coeff, func),...], 1:[...],</span>
<span class="sd">        ...} for all variables where i must be in [0..N-1] and tau &gt;= 0 with</span>
<span class="sd">        number of variables N. coeff must be a float and func a python</span>
<span class="sd">        callable of one argument.</span>
<span class="sd">    T : int</span>
<span class="sd">        Sample size.</span>
<span class="sd">    noises : list of callables or array, optional (default: &#39;np.random.randn&#39;)</span>
<span class="sd">        Random distribution function that is called with noises[j](T). If an array,</span>
<span class="sd">        it must be of shape ((transient_fraction + 1)*T, N).</span>
<span class="sd">    intervention : dict</span>
<span class="sd">        Dictionary of format: {1:np.array, ...} containing only keys of intervened</span>
<span class="sd">        variables with the value being the array of length T with interventional values.</span>
<span class="sd">        Set values to np.nan to leave specific time points of a variable un-intervened.</span>
<span class="sd">    intervention_type : str or dict</span>
<span class="sd">        Dictionary of format: {1:&#39;hard&#39;,  3:&#39;soft&#39;, ...} to specify whether intervention is </span>
<span class="sd">        hard (set value) or soft (add value) for variable j. If str, all interventions have </span>
<span class="sd">        the same type.</span>
<span class="sd">    transient_fraction : float</span>
<span class="sd">        Added percentage of T used as a transient. In total a realization of length</span>
<span class="sd">        (transient_fraction + 1)*T will be generated, but then transient_fraction*T will be</span>
<span class="sd">        cut off.</span>
<span class="sd">    seed : int, optional (default: None)</span>
<span class="sd">        Random seed.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    data : array-like</span>
<span class="sd">        Data generated from this process, shape (T, N).</span>
<span class="sd">    nonvalid : bool</span>
<span class="sd">        Indicates whether data has NaNs or infinities.</span>

<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">random_state</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">RandomState</span><span class="p">(</span><span class="n">seed</span><span class="p">)</span>

    <span class="n">N</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">links</span><span class="o">.</span><span class="n">keys</span><span class="p">())</span>
    <span class="k">if</span> <span class="n">noises</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
        <span class="n">noises</span> <span class="o">=</span> <span class="p">[</span><span class="n">random_state</span><span class="o">.</span><span class="n">randn</span> <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">)]</span>

    <span class="k">if</span> <span class="n">N</span> <span class="o">!=</span> <span class="nb">max</span><span class="p">(</span><span class="n">links</span><span class="o">.</span><span class="n">keys</span><span class="p">())</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;links keys must match N.&quot;</span><span class="p">)</span>

    <span class="k">if</span> <span class="nb">isinstance</span><span class="p">(</span><span class="n">noises</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">ndarray</span><span class="p">):</span>
        <span class="k">if</span> <span class="n">noises</span><span class="o">.</span><span class="n">shape</span> <span class="o">!=</span> <span class="p">(</span><span class="n">T</span> <span class="o">+</span> <span class="nb">int</span><span class="p">(</span><span class="n">math</span><span class="o">.</span><span class="n">floor</span><span class="p">(</span><span class="n">transient_fraction</span><span class="o">*</span><span class="n">T</span><span class="p">)),</span> <span class="n">N</span><span class="p">):</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;noises.shape must match ((transient_fraction + 1)*T, N).&quot;</span><span class="p">)</span>            
    <span class="k">else</span><span class="p">:</span>
        <span class="k">if</span> <span class="n">N</span> <span class="o">!=</span> <span class="nb">len</span><span class="p">(</span><span class="n">noises</span><span class="p">):</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;noises keys must match N.&quot;</span><span class="p">)</span>

    <span class="c1"># Check parameters</span>
    <span class="n">max_lag</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">contemp_dag</span> <span class="o">=</span> <span class="n">_Graph</span><span class="p">(</span><span class="n">N</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">link_props</span> <span class="ow">in</span> <span class="n">links</span><span class="p">[</span><span class="n">j</span><span class="p">]:</span>
            <span class="n">var</span><span class="p">,</span> <span class="n">lag</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
            <span class="n">coeff</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
            <span class="n">func</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
            <span class="k">if</span> <span class="n">lag</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span> <span class="n">contemp</span> <span class="o">=</span> <span class="kc">True</span>
            <span class="k">if</span> <span class="n">var</span> <span class="ow">not</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">):</span>
                <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;var must be in 0..</span><span class="si">{}</span><span class="s2">.&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">N</span><span class="o">-</span><span class="mi">1</span><span class="p">))</span>
            <span class="k">if</span> <span class="s1">&#39;float&#39;</span> <span class="ow">not</span> <span class="ow">in</span> <span class="nb">str</span><span class="p">(</span><span class="nb">type</span><span class="p">(</span><span class="n">coeff</span><span class="p">)):</span>
                <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;coeff must be float.&quot;</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">lag</span> <span class="o">&gt;</span> <span class="mi">0</span> <span class="ow">or</span> <span class="nb">type</span><span class="p">(</span><span class="n">lag</span><span class="p">)</span> <span class="o">!=</span> <span class="nb">int</span><span class="p">:</span>
                <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;lag must be non-positive int.&quot;</span><span class="p">)</span>
            <span class="n">max_lag</span> <span class="o">=</span> <span class="nb">max</span><span class="p">(</span><span class="n">max_lag</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">))</span>

            <span class="c1"># Create contemp DAG</span>
            <span class="k">if</span> <span class="n">var</span> <span class="o">!=</span> <span class="n">j</span> <span class="ow">and</span> <span class="n">lag</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
                <span class="n">contemp_dag</span><span class="o">.</span><span class="n">addEdge</span><span class="p">(</span><span class="n">var</span><span class="p">,</span> <span class="n">j</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">contemp_dag</span><span class="o">.</span><span class="n">isCyclic</span><span class="p">()</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span> 
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;Contemporaneous links must not contain cycle.&quot;</span><span class="p">)</span>

    <span class="n">causal_order</span> <span class="o">=</span> <span class="n">contemp_dag</span><span class="o">.</span><span class="n">topologicalSort</span><span class="p">()</span> 

    <span class="k">if</span> <span class="n">intervention</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
        <span class="k">if</span> <span class="n">intervention_type</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
            <span class="n">intervention_type</span> <span class="o">=</span> <span class="p">{</span><span class="n">j</span><span class="p">:</span><span class="s1">&#39;hard&#39;</span> <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">intervention</span><span class="p">}</span>
        <span class="k">elif</span> <span class="nb">isinstance</span><span class="p">(</span><span class="n">intervention_type</span><span class="p">,</span> <span class="nb">str</span><span class="p">):</span>
            <span class="n">intervention_type</span> <span class="o">=</span> <span class="p">{</span><span class="n">j</span><span class="p">:</span><span class="n">intervention_type</span> <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">intervention</span><span class="p">}</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">intervention</span><span class="o">.</span><span class="n">keys</span><span class="p">():</span>
            <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">intervention</span><span class="p">[</span><span class="n">j</span><span class="p">])</span> <span class="o">!=</span> <span class="n">T</span><span class="p">:</span>
                <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;intervention array for j=</span><span class="si">%s</span><span class="s2"> must be of length T = </span><span class="si">%d</span><span class="s2">&quot;</span> <span class="o">%</span><span class="p">(</span><span class="n">j</span><span class="p">,</span> <span class="n">T</span><span class="p">))</span>
            <span class="k">if</span> <span class="n">j</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">intervention_type</span><span class="o">.</span><span class="n">keys</span><span class="p">():</span>        
                <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;intervention_type dictionary must contain entry for </span><span class="si">%s</span><span class="s2">&quot;</span> <span class="o">%</span><span class="p">(</span><span class="n">j</span><span class="p">))</span>

    <span class="n">transient</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">math</span><span class="o">.</span><span class="n">floor</span><span class="p">(</span><span class="n">transient_fraction</span><span class="o">*</span><span class="n">T</span><span class="p">))</span>

    <span class="n">data</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">T</span><span class="o">+</span><span class="n">transient</span><span class="p">,</span> <span class="n">N</span><span class="p">),</span> <span class="n">dtype</span><span class="o">=</span><span class="s1">&#39;float32&#39;</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">):</span>
        <span class="k">if</span> <span class="nb">isinstance</span><span class="p">(</span><span class="n">noises</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">ndarray</span><span class="p">):</span>
            <span class="n">data</span><span class="p">[:,</span> <span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">noises</span><span class="p">[:,</span> <span class="n">j</span><span class="p">]</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">data</span><span class="p">[:,</span> <span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">noises</span><span class="p">[</span><span class="n">j</span><span class="p">](</span><span class="n">T</span><span class="o">+</span><span class="n">transient</span><span class="p">)</span>

    <span class="k">for</span> <span class="n">t</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">max_lag</span><span class="p">,</span> <span class="n">T</span><span class="o">+</span><span class="n">transient</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">causal_order</span><span class="p">:</span>

            <span class="k">if</span> <span class="p">(</span><span class="n">intervention</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span> <span class="ow">and</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">intervention</span> <span class="ow">and</span> <span class="n">t</span> <span class="o">&gt;=</span> <span class="n">transient</span>
                <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">intervention</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="n">t</span> <span class="o">-</span> <span class="n">transient</span><span class="p">])</span> <span class="o">==</span> <span class="kc">False</span><span class="p">):</span>
                <span class="k">if</span> <span class="n">intervention_type</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">==</span> <span class="s1">&#39;hard&#39;</span><span class="p">:</span>
                    <span class="n">data</span><span class="p">[</span><span class="n">t</span><span class="p">,</span> <span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">intervention</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="n">t</span> <span class="o">-</span> <span class="n">transient</span><span class="p">]</span>
                    <span class="c1"># Move to next j and skip link_props-loop from parents below </span>
                    <span class="k">continue</span>
                <span class="k">else</span><span class="p">:</span>
                    <span class="n">data</span><span class="p">[</span><span class="n">t</span><span class="p">,</span> <span class="n">j</span><span class="p">]</span> <span class="o">+=</span> <span class="n">intervention</span><span class="p">[</span><span class="n">j</span><span class="p">][</span><span class="n">t</span> <span class="o">-</span> <span class="n">transient</span><span class="p">]</span>

            <span class="c1"># This loop is only entered if intervention_type != &#39;hard&#39;</span>
            <span class="k">for</span> <span class="n">link_props</span> <span class="ow">in</span> <span class="n">links</span><span class="p">[</span><span class="n">j</span><span class="p">]:</span>
                <span class="n">var</span><span class="p">,</span> <span class="n">lag</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">coeff</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
                <span class="n">func</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
                <span class="n">data</span><span class="p">[</span><span class="n">t</span><span class="p">,</span> <span class="n">j</span><span class="p">]</span> <span class="o">+=</span> <span class="n">coeff</span> <span class="o">*</span> <span class="n">func</span><span class="p">(</span><span class="n">data</span><span class="p">[</span><span class="n">t</span> <span class="o">+</span> <span class="n">lag</span><span class="p">,</span> <span class="n">var</span><span class="p">])</span>

    <span class="n">data</span> <span class="o">=</span> <span class="n">data</span><span class="p">[</span><span class="n">transient</span><span class="p">:]</span>

    <span class="n">nonvalid</span> <span class="o">=</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">any</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">data</span><span class="p">))</span> <span class="ow">or</span> <span class="n">np</span><span class="o">.</span><span class="n">any</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">isinf</span><span class="p">(</span><span class="n">data</span><span class="p">)))</span>

    <span class="k">return</span> <span class="n">data</span><span class="p">,</span> <span class="n">nonvalid</span></div>

<span class="k">def</span> <span class="nf">_get_minmax_lag</span><span class="p">(</span><span class="n">links</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;Helper function to retrieve tau_min and tau_max from links.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">N</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">links</span><span class="p">)</span>

    <span class="c1"># Get maximum time lag</span>
    <span class="n">min_lag</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">inf</span>
    <span class="n">max_lag</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">link_props</span> <span class="ow">in</span> <span class="n">links</span><span class="p">[</span><span class="n">j</span><span class="p">]:</span>
            <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">link_props</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">2</span><span class="p">:</span>
                <span class="n">var</span><span class="p">,</span> <span class="n">lag</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">coeff</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
                <span class="c1"># func = link_props[2]</span>
                <span class="k">if</span> <span class="ow">not</span> <span class="nb">isinstance</span><span class="p">(</span><span class="n">coeff</span><span class="p">,</span> <span class="nb">float</span><span class="p">)</span> <span class="ow">or</span> <span class="n">coeff</span> <span class="o">!=</span> <span class="mf">0.</span><span class="p">:</span>
                    <span class="n">min_lag</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="n">min_lag</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">))</span>
                    <span class="n">max_lag</span> <span class="o">=</span> <span class="nb">max</span><span class="p">(</span><span class="n">max_lag</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">))</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">var</span><span class="p">,</span> <span class="n">lag</span> <span class="o">=</span> <span class="n">link_props</span>
                <span class="n">min_lag</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="n">min_lag</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">))</span>
                <span class="n">max_lag</span> <span class="o">=</span> <span class="nb">max</span><span class="p">(</span><span class="n">max_lag</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">))</span>   

    <span class="k">return</span> <span class="n">min_lag</span><span class="p">,</span> <span class="n">max_lag</span>

<span class="k">def</span> <span class="nf">_get_parents</span><span class="p">(</span><span class="n">links</span><span class="p">,</span> <span class="n">exclude_contemp</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;Helper function to parents from links</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">N</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">links</span><span class="p">)</span>

    <span class="c1"># Get maximum time lag</span>
    <span class="n">parents</span> <span class="o">=</span> <span class="p">{}</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">):</span>
        <span class="n">parents</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="p">[]</span>
        <span class="k">for</span> <span class="n">link_props</span> <span class="ow">in</span> <span class="n">links</span><span class="p">[</span><span class="n">j</span><span class="p">]:</span>
            <span class="n">var</span><span class="p">,</span> <span class="n">lag</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
            <span class="n">coeff</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
            <span class="c1"># func = link_props[2]</span>
            <span class="k">if</span> <span class="n">coeff</span> <span class="o">!=</span> <span class="mf">0.</span><span class="p">:</span>
                <span class="k">if</span> <span class="ow">not</span> <span class="p">(</span><span class="n">exclude_contemp</span> <span class="ow">and</span> <span class="n">lag</span> <span class="o">==</span> <span class="mi">0</span><span class="p">):</span>
                    <span class="n">parents</span><span class="p">[</span><span class="n">j</span><span class="p">]</span><span class="o">.</span><span class="n">append</span><span class="p">((</span><span class="n">var</span><span class="p">,</span> <span class="n">lag</span><span class="p">))</span>

    <span class="k">return</span> <span class="n">parents</span>

<span class="k">def</span> <span class="nf">_get_children</span><span class="p">(</span><span class="n">parents</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;Helper function to children from parents</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">N</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">parents</span><span class="p">)</span>
    <span class="n">children</span> <span class="o">=</span> <span class="nb">dict</span><span class="p">([(</span><span class="n">j</span><span class="p">,</span> <span class="p">[])</span> <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">)])</span>

    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">par</span> <span class="ow">in</span> <span class="n">parents</span><span class="p">[</span><span class="n">j</span><span class="p">]:</span>
            <span class="n">i</span><span class="p">,</span> <span class="n">tau</span> <span class="o">=</span> <span class="n">par</span>
            <span class="n">children</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">append</span><span class="p">((</span><span class="n">j</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">tau</span><span class="p">)))</span>

    <span class="k">return</span> <span class="n">children</span>

<div class="viewcode-block" id="links_to_graph"><a class="viewcode-back" href="../../../index.html#tigramite.toymodels.structural_causal_processes.links_to_graph">[docs]</a><span class="k">def</span> <span class="nf">links_to_graph</span><span class="p">(</span><span class="n">links</span><span class="p">,</span> <span class="n">tau_max</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;Helper function to convert dictionary of links to graph array format.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ---------</span>
<span class="sd">    links : dict</span>
<span class="sd">        Dictionary of form {0:[((0, -1), coeff, func), ...], 1:[...], ...}.</span>
<span class="sd">        Also format {0:[(0, -1), ...], 1:[...], ...} is allowed.</span>
<span class="sd">    tau_max : int or None</span>
<span class="sd">        Maximum lag. If None, the maximum lag in links is used.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    graph : array of shape (N, N, tau_max+1)</span>
<span class="sd">        Matrix format of graph with 1 for true links and 0 else.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">N</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">links</span><span class="p">)</span>

    <span class="c1"># Get maximum time lag</span>
    <span class="n">min_lag</span><span class="p">,</span> <span class="n">max_lag</span> <span class="o">=</span> <span class="n">_get_minmax_lag</span><span class="p">(</span><span class="n">links</span><span class="p">)</span>

    <span class="c1"># Set maximum lag</span>
    <span class="k">if</span> <span class="n">tau_max</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
        <span class="n">tau_max</span> <span class="o">=</span> <span class="n">max_lag</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="k">if</span> <span class="n">max_lag</span> <span class="o">&gt;</span> <span class="n">tau_max</span><span class="p">:</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;tau_max is smaller than maximum lag = </span><span class="si">%d</span><span class="s2"> &quot;</span>
                             <span class="s2">&quot;found in links, use tau_max=None or larger &quot;</span>
                             <span class="s2">&quot;value&quot;</span> <span class="o">%</span> <span class="n">max_lag</span><span class="p">)</span>

    <span class="n">graph</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">N</span><span class="p">,</span> <span class="n">N</span><span class="p">,</span> <span class="n">tau_max</span> <span class="o">+</span> <span class="mi">1</span><span class="p">),</span> <span class="n">dtype</span><span class="o">=</span><span class="s1">&#39;&lt;U3&#39;</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">links</span><span class="o">.</span><span class="n">keys</span><span class="p">():</span>
        <span class="k">for</span> <span class="n">link_props</span> <span class="ow">in</span> <span class="n">links</span><span class="p">[</span><span class="n">j</span><span class="p">]:</span>
            <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">link_props</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">2</span><span class="p">:</span>
                <span class="n">var</span><span class="p">,</span> <span class="n">lag</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">coeff</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
                <span class="k">if</span> <span class="n">coeff</span> <span class="o">!=</span> <span class="mf">0.</span><span class="p">:</span>
                    <span class="n">graph</span><span class="p">[</span><span class="n">var</span><span class="p">,</span> <span class="n">j</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">)]</span> <span class="o">=</span> <span class="s2">&quot;--&gt;&quot;</span>
                    <span class="k">if</span> <span class="n">lag</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
                        <span class="n">graph</span><span class="p">[</span><span class="n">j</span><span class="p">,</span> <span class="n">var</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="s2">&quot;&lt;--&quot;</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">var</span><span class="p">,</span> <span class="n">lag</span> <span class="o">=</span> <span class="n">link_props</span>
                <span class="n">graph</span><span class="p">[</span><span class="n">var</span><span class="p">,</span> <span class="n">j</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">)]</span> <span class="o">=</span> <span class="s2">&quot;--&gt;&quot;</span>
                <span class="k">if</span> <span class="n">lag</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
                    <span class="n">graph</span><span class="p">[</span><span class="n">j</span><span class="p">,</span> <span class="n">var</span><span class="p">,</span> <span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="s2">&quot;&lt;--&quot;</span>

    <span class="k">return</span> <span class="n">graph</span></div>

<div class="viewcode-block" id="dag_to_links"><a class="viewcode-back" href="../../../index.html#tigramite.toymodels.structural_causal_processes.dag_to_links">[docs]</a><span class="k">def</span> <span class="nf">dag_to_links</span><span class="p">(</span><span class="n">dag</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;Helper function to convert DAG graph to dictionary of parents.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ---------</span>
<span class="sd">    dag : array of shape (N, N, tau_max+1)</span>
<span class="sd">        Matrix format of graph in string format. Must be DAG.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    parents : dict</span>
<span class="sd">        Dictionary of form {0:[(0, -1), ...], 1:[...], ...}.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">N</span> <span class="o">=</span> <span class="n">dag</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>

    <span class="n">parents</span> <span class="o">=</span> <span class="nb">dict</span><span class="p">([(</span><span class="n">j</span><span class="p">,</span> <span class="p">[])</span> <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">)])</span>

    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">any</span><span class="p">(</span><span class="n">dag</span><span class="o">==</span><span class="s1">&#39;o-o&#39;</span><span class="p">)</span> <span class="ow">or</span> <span class="n">np</span><span class="o">.</span><span class="n">any</span><span class="p">(</span><span class="n">dag</span><span class="o">==</span><span class="s1">&#39;x-x&#39;</span><span class="p">):</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;graph must be DAG.&quot;</span><span class="p">)</span>

    <span class="k">for</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">j</span><span class="p">,</span> <span class="n">tau</span><span class="p">)</span> <span class="ow">in</span> <span class="nb">zip</span><span class="p">(</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">where</span><span class="p">(</span><span class="n">dag</span><span class="o">==</span><span class="s1">&#39;--&gt;&#39;</span><span class="p">)):</span>
        <span class="n">parents</span><span class="p">[</span><span class="n">j</span><span class="p">]</span><span class="o">.</span><span class="n">append</span><span class="p">((</span><span class="n">i</span><span class="p">,</span> <span class="o">-</span><span class="n">tau</span><span class="p">))</span>

    <span class="k">return</span> <span class="n">parents</span></div>

<div class="viewcode-block" id="generate_structural_causal_process"><a class="viewcode-back" href="../../../index.html#tigramite.toymodels.structural_causal_processes.generate_structural_causal_process">[docs]</a><span class="k">def</span> <span class="nf">generate_structural_causal_process</span><span class="p">(</span>
                <span class="n">N</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> 
                <span class="n">L</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> 
                <span class="n">dependency_funcs</span><span class="o">=</span><span class="p">[</span><span class="s1">&#39;linear&#39;</span><span class="p">],</span> 
                <span class="n">dependency_coeffs</span><span class="o">=</span><span class="p">[</span><span class="o">-</span><span class="mf">0.5</span><span class="p">,</span> <span class="mf">0.5</span><span class="p">],</span> 
                <span class="n">auto_coeffs</span><span class="o">=</span><span class="p">[</span><span class="mf">0.5</span><span class="p">,</span> <span class="mf">0.7</span><span class="p">],</span> 
                <span class="n">contemp_fraction</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span>
                <span class="n">max_lag</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> 
                <span class="n">noise_dists</span><span class="o">=</span><span class="p">[</span><span class="s1">&#39;gaussian&#39;</span><span class="p">],</span>
                <span class="n">noise_means</span><span class="o">=</span><span class="p">[</span><span class="mf">0.</span><span class="p">],</span>
                <span class="n">noise_sigmas</span><span class="o">=</span><span class="p">[</span><span class="mf">0.5</span><span class="p">,</span> <span class="mf">2.</span><span class="p">],</span>
                <span class="n">noise_seed</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span>
                <span class="n">seed</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;&quot;Randomly generates a structural causal process based on input characteristics.</span>

<span class="sd">    The process has the form </span>

<span class="sd">    .. math:: X^j_t = \\eta^j_t + a^j X^j_{t-1} + \\sum_{X^i_{t-\\tau}\\in pa(X^j_t)}</span>
<span class="sd">              c^i_{\\tau} f^i_{\\tau}(X^i_{t-\\tau})</span>

<span class="sd">    where ``j = 1, ..., N``. Here the properties of :math:`\\eta^j_t` are</span>
<span class="sd">    randomly frawn from the noise parameters (see below), :math:`pa</span>
<span class="sd">    (X^j_t)` are the causal parents drawn randomly such that in total ``L``</span>
<span class="sd">    links occur out of which ``contemp_fraction`` are contemporaneous and</span>
<span class="sd">    their time lags are drawn from ``[0 or 1..max_lag]``, the</span>
<span class="sd">    coefficients :math:`c^i_{\\tau}` are drawn from</span>
<span class="sd">    ``dependency_coeffs``, :math:`a^j` are drawn from ``auto_coeffs``,</span>
<span class="sd">    and :math:`f^i_{\\tau}` are drawn from ``dependency_funcs``.</span>

<span class="sd">    The returned dictionary links has the format </span>
<span class="sd">    ``{0:[((i, -tau), coeff, func),...], 1:[...], ...}`` </span>
<span class="sd">    where ``func`` can be an arbitrary (nonlinear) function provided</span>
<span class="sd">    as a python callable with one argument and coeff is the multiplication</span>
<span class="sd">    factor. The noise distributions of :math:`\\eta^j` are returned in</span>
<span class="sd">    ``noises``, see specifics below.</span>

<span class="sd">    The process might be non-stationary. In case of asymptotically linear</span>
<span class="sd">    dependency functions and no contemporaneous links this can be checked with</span>
<span class="sd">    ``check_stationarity(...)``. Otherwise check by generating a large sample</span>
<span class="sd">    and test for np.inf.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ---------</span>
<span class="sd">    N : int</span>
<span class="sd">        Number of variables.</span>
<span class="sd">    L : int</span>
<span class="sd">        Number of cross-links between two different variables.</span>
<span class="sd">    dependency_funcs : list</span>
<span class="sd">        List of callables or strings &#39;linear&#39; or &#39;nonlinear&#39; for a linear and a specific nonlinear function</span>
<span class="sd">        that is asymptotically linear.</span>
<span class="sd">    dependency_coeffs : list</span>
<span class="sd">        List of floats from which the coupling coefficients are randomly drawn.</span>
<span class="sd">    auto_coeffs : list</span>
<span class="sd">        List of floats from which the lag-1 autodependencies are randomly drawn.</span>
<span class="sd">    contemp_fraction : float [0., 1]</span>
<span class="sd">        Fraction of the L links that are contemporaneous (lag zero).</span>
<span class="sd">    max_lag : int</span>
<span class="sd">        Maximum lag from which the time lags of links are drawn.</span>
<span class="sd">    noise_dists : list</span>
<span class="sd">        List of noise functions. Either in</span>
<span class="sd">        {&#39;gaussian&#39;, &#39;weibull&#39;, &#39;uniform&#39;} or user-specified, in which case</span>
<span class="sd">        it must be parametrized just by the size parameter. E.g. def beta</span>
<span class="sd">        (T): return np.random.beta(a=1, b=0.5, T)</span>
<span class="sd">    noise_means : list</span>
<span class="sd">        Noise mean. Only used for noise in {&#39;gaussian&#39;, &#39;weibull&#39;, &#39;uniform&#39;}.</span>
<span class="sd">    noise_sigmas : list</span>
<span class="sd">        Noise standard deviation. Only used for noise in {&#39;gaussian&#39;, &#39;weibull&#39;, &#39;uniform&#39;}.   </span>
<span class="sd">    seed : int</span>
<span class="sd">        Random seed to draw the above random functions from.</span>
<span class="sd">    noise_seed : int</span>
<span class="sd">        Random seed for noise function random generator.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    links : dict</span>
<span class="sd">        Dictionary of form {0:[((0, -1), coeff, func), ...], 1:[...], ...}.</span>
<span class="sd">    noises : list</span>
<span class="sd">        List of N noise functions to call by noise(T) where T is the time series length.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    
    <span class="c1"># Init random states</span>
    <span class="n">random_state</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">RandomState</span><span class="p">(</span><span class="n">seed</span><span class="p">)</span>
    <span class="n">random_state_noise</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">RandomState</span><span class="p">(</span><span class="n">noise_seed</span><span class="p">)</span>

    <span class="k">def</span> <span class="nf">linear</span><span class="p">(</span><span class="n">x</span><span class="p">):</span> <span class="k">return</span> <span class="n">x</span>
    <span class="k">def</span> <span class="nf">nonlinear</span><span class="p">(</span><span class="n">x</span><span class="p">):</span> <span class="k">return</span> <span class="p">(</span><span class="n">x</span> <span class="o">+</span> <span class="mf">5.</span> <span class="o">*</span> <span class="n">x</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="mf">20.</span><span class="p">))</span>

    <span class="k">if</span> <span class="n">max_lag</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">contemp_fraction</span> <span class="o">=</span> <span class="mf">1.</span>

    <span class="k">if</span> <span class="n">contemp_fraction</span> <span class="o">&gt;</span> <span class="mf">0.</span><span class="p">:</span>
        <span class="n">ordered_pairs</span> <span class="o">=</span> <span class="nb">list</span><span class="p">(</span><span class="n">itertools</span><span class="o">.</span><span class="n">combinations</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">),</span> <span class="mi">2</span><span class="p">))</span>
        <span class="n">max_poss_links</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="n">L</span><span class="p">,</span> <span class="nb">len</span><span class="p">(</span><span class="n">ordered_pairs</span><span class="p">))</span>
        <span class="n">L_contemp</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">contemp_fraction</span><span class="o">*</span><span class="n">max_poss_links</span><span class="p">)</span>
        <span class="n">L_lagged</span> <span class="o">=</span> <span class="n">max_poss_links</span> <span class="o">-</span> <span class="n">L_contemp</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">L_lagged</span> <span class="o">=</span> <span class="n">L</span>
        <span class="n">L_contemp</span> <span class="o">=</span> <span class="mi">0</span>

    <span class="c1"># Random order</span>
    <span class="n">causal_order</span> <span class="o">=</span> <span class="nb">list</span><span class="p">(</span><span class="n">random_state</span><span class="o">.</span><span class="n">permutation</span><span class="p">(</span><span class="n">N</span><span class="p">))</span>

    <span class="c1"># Init link dict</span>
    <span class="n">links</span> <span class="o">=</span> <span class="nb">dict</span><span class="p">([(</span><span class="n">i</span><span class="p">,</span> <span class="p">[])</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">)])</span>

    <span class="c1"># Generate auto-dependencies at lag 1</span>
    <span class="k">if</span> <span class="n">max_lag</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">causal_order</span><span class="p">:</span>
            <span class="n">a</span> <span class="o">=</span> <span class="n">random_state</span><span class="o">.</span><span class="n">choice</span><span class="p">(</span><span class="n">auto_coeffs</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">a</span> <span class="o">!=</span> <span class="mf">0.</span><span class="p">:</span>
                <span class="n">links</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">append</span><span class="p">(((</span><span class="nb">int</span><span class="p">(</span><span class="n">i</span><span class="p">),</span> <span class="o">-</span><span class="mi">1</span><span class="p">),</span> <span class="nb">float</span><span class="p">(</span><span class="n">a</span><span class="p">),</span> <span class="n">linear</span><span class="p">))</span>

    <span class="c1"># Non-cyclic contemp random pairs of links such that</span>
    <span class="c1"># index of cause &lt; index of effect</span>
    <span class="c1"># Take up to (!) L_contemp links</span>
    <span class="n">ordered_pairs</span> <span class="o">=</span> <span class="nb">list</span><span class="p">(</span><span class="n">itertools</span><span class="o">.</span><span class="n">combinations</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">),</span> <span class="mi">2</span><span class="p">))</span>
    <span class="n">random_state</span><span class="o">.</span><span class="n">shuffle</span><span class="p">(</span><span class="n">ordered_pairs</span><span class="p">)</span>
    <span class="n">contemp_links</span> <span class="o">=</span> <span class="p">[(</span><span class="n">causal_order</span><span class="p">[</span><span class="n">pair</span><span class="p">[</span><span class="mi">0</span><span class="p">]],</span> <span class="n">causal_order</span><span class="p">[</span><span class="n">pair</span><span class="p">[</span><span class="mi">1</span><span class="p">]])</span> 
                    <span class="k">for</span> <span class="n">pair</span> <span class="ow">in</span> <span class="n">ordered_pairs</span><span class="p">[:</span><span class="n">L_contemp</span><span class="p">]]</span>

    <span class="c1"># Possibly cyclic lagged random pairs of links </span>
    <span class="c1"># where we remove already chosen contemp links</span>
    <span class="c1"># Take up to (!) L_contemp links</span>
    <span class="n">unordered_pairs</span> <span class="o">=</span> <span class="nb">list</span><span class="p">(</span><span class="n">itertools</span><span class="o">.</span><span class="n">permutations</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">),</span> <span class="mi">2</span><span class="p">))</span>
    <span class="n">unordered_pairs</span> <span class="o">=</span> <span class="nb">list</span><span class="p">(</span><span class="nb">set</span><span class="p">(</span><span class="n">unordered_pairs</span><span class="p">)</span> <span class="o">-</span> <span class="nb">set</span><span class="p">(</span><span class="n">ordered_pairs</span><span class="p">[:</span><span class="n">L_contemp</span><span class="p">]))</span>
    <span class="n">random_state</span><span class="o">.</span><span class="n">shuffle</span><span class="p">(</span><span class="n">unordered_pairs</span><span class="p">)</span>
    <span class="n">lagged_links</span> <span class="o">=</span> <span class="p">[(</span><span class="n">causal_order</span><span class="p">[</span><span class="n">pair</span><span class="p">[</span><span class="mi">0</span><span class="p">]],</span> <span class="n">causal_order</span><span class="p">[</span><span class="n">pair</span><span class="p">[</span><span class="mi">1</span><span class="p">]])</span> 
                    <span class="k">for</span> <span class="n">pair</span> <span class="ow">in</span> <span class="n">unordered_pairs</span><span class="p">[:</span><span class="n">L_lagged</span><span class="p">]]</span>

    <span class="n">chosen_links</span> <span class="o">=</span> <span class="n">lagged_links</span> <span class="o">+</span> <span class="n">contemp_links</span>

    <span class="c1"># Populate links</span>
    <span class="k">for</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">j</span><span class="p">)</span> <span class="ow">in</span> <span class="n">chosen_links</span><span class="p">:</span>

        <span class="c1"># Choose lag</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">j</span><span class="p">)</span> <span class="ow">in</span> <span class="n">contemp_links</span><span class="p">:</span>
            <span class="n">tau</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">tau</span> <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">random_state</span><span class="o">.</span><span class="n">randint</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="n">max_lag</span><span class="o">+</span><span class="mi">1</span><span class="p">))</span>

        <span class="c1"># Choose dependency</span>
        <span class="n">c</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">random_state</span><span class="o">.</span><span class="n">choice</span><span class="p">(</span><span class="n">dependency_coeffs</span><span class="p">))</span>
        <span class="k">if</span> <span class="n">c</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">func</span> <span class="o">=</span> <span class="n">random_state</span><span class="o">.</span><span class="n">choice</span><span class="p">(</span><span class="n">dependency_funcs</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">func</span> <span class="o">==</span> <span class="s1">&#39;linear&#39;</span><span class="p">:</span>
                <span class="n">func</span> <span class="o">=</span> <span class="n">linear</span>
            <span class="k">elif</span> <span class="n">func</span> <span class="o">==</span> <span class="s1">&#39;nonlinear&#39;</span><span class="p">:</span>
                <span class="n">func</span> <span class="o">=</span> <span class="n">nonlinear</span>

            <span class="n">links</span><span class="p">[</span><span class="n">j</span><span class="p">]</span><span class="o">.</span><span class="n">append</span><span class="p">(((</span><span class="nb">int</span><span class="p">(</span><span class="n">i</span><span class="p">),</span> <span class="o">-</span><span class="n">tau</span><span class="p">),</span> <span class="n">c</span><span class="p">,</span> <span class="n">func</span><span class="p">))</span>

    <span class="c1"># Now generate noise functions</span>
    <span class="c1"># Either choose among pre-defined noise types or supply your own</span>
    <span class="k">class</span> <span class="nc">NoiseModel</span><span class="p">:</span>
        <span class="k">def</span> <span class="fm">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">mean</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span> <span class="n">sigma</span><span class="o">=</span><span class="mf">1.</span><span class="p">):</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">mean</span> <span class="o">=</span> <span class="n">mean</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">sigma</span> <span class="o">=</span> <span class="n">sigma</span>
        <span class="k">def</span> <span class="nf">gaussian</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">T</span><span class="p">):</span>
            <span class="c1"># Get zero-mean unit variance gaussian distribution</span>
            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">mean</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">sigma</span><span class="o">*</span><span class="n">random_state_noise</span><span class="o">.</span><span class="n">randn</span><span class="p">(</span><span class="n">T</span><span class="p">)</span>
        <span class="k">def</span> <span class="nf">weibull</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">T</span><span class="p">):</span> 
            <span class="c1"># Get zero-mean sigma variance weibull distribution</span>
            <span class="n">a</span> <span class="o">=</span> <span class="mi">2</span>
            <span class="n">mean</span> <span class="o">=</span> <span class="n">scipy</span><span class="o">.</span><span class="n">special</span><span class="o">.</span><span class="n">gamma</span><span class="p">(</span><span class="mf">1.</span><span class="o">/</span><span class="n">a</span> <span class="o">+</span> <span class="mi">1</span><span class="p">)</span>
            <span class="n">variance</span> <span class="o">=</span> <span class="n">scipy</span><span class="o">.</span><span class="n">special</span><span class="o">.</span><span class="n">gamma</span><span class="p">(</span><span class="mf">2.</span><span class="o">/</span><span class="n">a</span> <span class="o">+</span> <span class="mi">1</span><span class="p">)</span> <span class="o">-</span> <span class="n">scipy</span><span class="o">.</span><span class="n">special</span><span class="o">.</span><span class="n">gamma</span><span class="p">(</span><span class="mf">1.</span><span class="o">/</span><span class="n">a</span> <span class="o">+</span> <span class="mi">1</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span>
            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">mean</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">sigma</span><span class="o">*</span><span class="p">(</span><span class="n">random_state_noise</span><span class="o">.</span><span class="n">weibull</span><span class="p">(</span><span class="n">a</span><span class="o">=</span><span class="n">a</span><span class="p">,</span> <span class="n">size</span><span class="o">=</span><span class="n">T</span><span class="p">)</span> <span class="o">-</span> <span class="n">mean</span><span class="p">)</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">variance</span><span class="p">)</span>
        <span class="k">def</span> <span class="nf">uniform</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">T</span><span class="p">):</span> 
            <span class="c1"># Get zero-mean sigma variance uniform distribution</span>
            <span class="n">mean</span> <span class="o">=</span> <span class="mf">0.5</span>
            <span class="n">variance</span> <span class="o">=</span> <span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span>
            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">mean</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">sigma</span><span class="o">*</span><span class="p">(</span><span class="n">random_state_noise</span><span class="o">.</span><span class="n">uniform</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">T</span><span class="p">)</span> <span class="o">-</span> <span class="n">mean</span><span class="p">)</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">variance</span><span class="p">)</span>

    <span class="n">noises</span> <span class="o">=</span> <span class="p">[]</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">links</span><span class="p">:</span>
        <span class="n">noise_dist</span> <span class="o">=</span> <span class="n">random_state</span><span class="o">.</span><span class="n">choice</span><span class="p">(</span><span class="n">noise_dists</span><span class="p">)</span>
        <span class="n">noise_mean</span> <span class="o">=</span> <span class="n">random_state</span><span class="o">.</span><span class="n">choice</span><span class="p">(</span><span class="n">noise_means</span><span class="p">)</span>
        <span class="n">noise_sigma</span> <span class="o">=</span> <span class="n">random_state</span><span class="o">.</span><span class="n">choice</span><span class="p">(</span><span class="n">noise_sigmas</span><span class="p">)</span>

        <span class="k">if</span> <span class="n">noise_dist</span> <span class="ow">in</span> <span class="p">[</span><span class="s1">&#39;gaussian&#39;</span><span class="p">,</span> <span class="s1">&#39;weibull&#39;</span><span class="p">,</span> <span class="s1">&#39;uniform&#39;</span><span class="p">]:</span>
            <span class="n">noise</span> <span class="o">=</span> <span class="nb">getattr</span><span class="p">(</span><span class="n">NoiseModel</span><span class="p">(</span><span class="n">mean</span> <span class="o">=</span> <span class="n">noise_mean</span><span class="p">,</span> <span class="n">sigma</span> <span class="o">=</span> <span class="n">noise_sigma</span><span class="p">),</span> <span class="n">noise_dist</span><span class="p">)</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">noise</span> <span class="o">=</span> <span class="n">noise_dist</span>
        
        <span class="n">noises</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">noise</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">links</span><span class="p">,</span> <span class="n">noises</span></div>

<div class="viewcode-block" id="check_stationarity"><a class="viewcode-back" href="../../../index.html#tigramite.toymodels.structural_causal_processes.check_stationarity">[docs]</a><span class="k">def</span> <span class="nf">check_stationarity</span><span class="p">(</span><span class="n">links</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;Returns stationarity according to a unit root test.</span>

<span class="sd">    Assumes an at least asymptotically linear vector autoregressive process</span>
<span class="sd">    without contemporaneous links.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ---------</span>
<span class="sd">    links : dict</span>
<span class="sd">        Dictionary of form {0:[((0, -1), coeff, func), ...], 1:[...], ...}.</span>
<span class="sd">        Also format {0:[(0, -1), ...], 1:[...], ...} is allowed.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    stationary : bool</span>
<span class="sd">        True if VAR process is stationary.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">N</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">links</span><span class="p">)</span>
    <span class="c1"># Check parameters</span>
    <span class="n">max_lag</span> <span class="o">=</span> <span class="mi">0</span>

    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">link_props</span> <span class="ow">in</span> <span class="n">links</span><span class="p">[</span><span class="n">j</span><span class="p">]:</span>
            <span class="n">var</span><span class="p">,</span> <span class="n">lag</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
            <span class="c1"># coeff = link_props[1]</span>
            <span class="c1"># coupling = link_props[2]</span>

            <span class="n">max_lag</span> <span class="o">=</span> <span class="nb">max</span><span class="p">(</span><span class="n">max_lag</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">))</span>

    <span class="n">graph</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">N</span><span class="p">,</span><span class="n">N</span><span class="p">,</span><span class="n">max_lag</span><span class="p">))</span>
    <span class="n">couplings</span> <span class="o">=</span> <span class="p">[]</span>

    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">N</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">link_props</span> <span class="ow">in</span> <span class="n">links</span><span class="p">[</span><span class="n">j</span><span class="p">]:</span>
            <span class="n">var</span><span class="p">,</span> <span class="n">lag</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
            <span class="n">coeff</span>    <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
            <span class="n">coupling</span> <span class="o">=</span> <span class="n">link_props</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
            <span class="k">if</span> <span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">)</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
                <span class="n">graph</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">var</span><span class="p">,</span><span class="nb">abs</span><span class="p">(</span><span class="n">lag</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">coeff</span>
            <span class="n">couplings</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">coupling</span><span class="p">)</span>

    <span class="n">stabmat</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">N</span><span class="o">*</span><span class="n">max_lag</span><span class="p">,</span><span class="n">N</span><span class="o">*</span><span class="n">max_lag</span><span class="p">))</span>
    <span class="n">index</span> <span class="o">=</span> <span class="mi">0</span>

    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="n">N</span><span class="o">*</span><span class="n">max_lag</span><span class="p">,</span><span class="n">N</span><span class="p">):</span>
        <span class="n">stabmat</span><span class="p">[:</span><span class="n">N</span><span class="p">,</span><span class="n">i</span><span class="p">:</span><span class="n">i</span><span class="o">+</span><span class="n">N</span><span class="p">]</span> <span class="o">=</span> <span class="n">graph</span><span class="p">[:,:,</span><span class="n">index</span><span class="p">]</span>
        <span class="k">if</span> <span class="n">index</span> <span class="o">&lt;</span> <span class="n">max_lag</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span>
            <span class="n">stabmat</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="n">N</span><span class="p">:</span><span class="n">i</span><span class="o">+</span><span class="mi">2</span><span class="o">*</span><span class="n">N</span><span class="p">,</span><span class="n">i</span><span class="p">:</span><span class="n">i</span><span class="o">+</span><span class="n">N</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">identity</span><span class="p">(</span><span class="n">N</span><span class="p">)</span>
        <span class="n">index</span> <span class="o">+=</span> <span class="mi">1</span>

    <span class="n">eig</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">eig</span><span class="p">(</span><span class="n">stabmat</span><span class="p">)[</span><span class="mi">0</span><span class="p">]</span>

    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">all</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">eig</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1.</span><span class="p">):</span>
        <span class="n">stationary</span> <span class="o">=</span> <span class="kc">True</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">stationary</span> <span class="o">=</span> <span class="kc">False</span>

    <span class="k">return</span> <span class="n">stationary</span></div>
    <span class="c1"># if len(eig) == 0:</span>
    <span class="c1">#     return stationary, 0.</span>
    <span class="c1"># else:</span>
    <span class="c1">#     return stationary, np.abs(eig).max()</span>

<span class="k">class</span> <span class="nc">_Logger</span><span class="p">(</span><span class="nb">object</span><span class="p">):</span>
<span class="w">    </span><span class="sd">&quot;&quot;&quot;Class to append print output to a string which can be saved&quot;&quot;&quot;</span>
    <span class="k">def</span> <span class="fm">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">terminal</span> <span class="o">=</span> <span class="n">sys</span><span class="o">.</span><span class="n">stdout</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">log</span> <span class="o">=</span> <span class="s2">&quot;&quot;</span>       <span class="c1"># open(&quot;log.dat&quot;, &quot;a&quot;)</span>

    <span class="k">def</span> <span class="nf">write</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">message</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">terminal</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="n">message</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">log</span> <span class="o">+=</span> <span class="n">message</span>  <span class="c1"># .write(message)</span>


<span class="k">if</span> <span class="vm">__name__</span> <span class="o">==</span> <span class="s1">&#39;__main__&#39;</span><span class="p">:</span>
    
    <span class="c1">## Generate some time series from a structural causal process</span>
    <span class="k">def</span> <span class="nf">lin_f</span><span class="p">(</span><span class="n">x</span><span class="p">):</span> <span class="k">return</span> <span class="n">x</span>
    <span class="k">def</span> <span class="nf">nonlin_f</span><span class="p">(</span><span class="n">x</span><span class="p">):</span> <span class="k">return</span> <span class="p">(</span><span class="n">x</span> <span class="o">+</span> <span class="mf">5.</span> <span class="o">*</span> <span class="n">x</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="mf">20.</span><span class="p">))</span>

    <span class="n">links</span><span class="p">,</span> <span class="n">noises</span> <span class="o">=</span> <span class="n">generate_structural_causal_process</span><span class="p">()</span>
    
    <span class="n">data</span><span class="p">,</span> <span class="n">nonstat</span> <span class="o">=</span> <span class="n">structural_causal_process</span><span class="p">(</span><span class="n">links</span><span class="p">,</span>
             <span class="n">T</span><span class="o">=</span><span class="mi">100</span><span class="p">,</span> <span class="n">noises</span><span class="o">=</span><span class="n">noises</span><span class="p">)</span>
    <span class="nb">print</span><span class="p">(</span><span class="n">data</span><span class="p">)</span>

    <span class="c1"># links = {0: [((0, -1), 0.9, lin_f)],</span>
    <span class="c1">#          1: [((1, -1), 0.8, lin_f), ((0, -1), 0.3, nonlin_f)],</span>
    <span class="c1">#          2: [((2, -1), 0.7, lin_f), ((1, 0), -0.2, lin_f)],</span>
    <span class="c1">#          }</span>
    <span class="c1"># noises = [np.random.randn, np.random.randn, np.random.randn]</span>

    <span class="c1"># data, nonstat = structural_causal_process(links,</span>
    <span class="c1">#  T=100, noises=noises)</span>

</pre></div>

          </div>
          
        </div>
      </div>
      <div class="sphinxsidebar" role="navigation" aria-label="main navigation">
        <div class="sphinxsidebarwrapper">
<h1 class="logo"><a href="../../../index.html">Tigramite</a></h1>








<h3>Navigation</h3>

<div class="relations">
<h3>Related Topics</h3>
<ul>
  <li><a href="../../../index.html">Documentation overview</a><ul>
  <li><a href="../../index.html">Module code</a><ul>
  </ul></li>
  </ul></li>
</ul>
</div>
<div id="searchbox" style="display: none" role="search">
  <h3 id="searchlabel">Quick search</h3>
    <div class="searchformwrapper">
    <form class="search" action="../../../search.html" method="get">
      <input type="text" name="q" aria-labelledby="searchlabel" autocomplete="off" autocorrect="off" autocapitalize="off" spellcheck="false"/>
      <input type="submit" value="Go" />
    </form>
    </div>
</div>
<script>document.getElementById('searchbox').style.display = "block"</script>








        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="footer">
      &copy;2023, Jakob Runge.
      
      |
      Powered by <a href="http://sphinx-doc.org/">Sphinx 5.0.2</a>
      &amp; <a href="https://github.com/bitprophet/alabaster">Alabaster 0.7.12</a>
      
    </div>

    

    
  </body>
</html>